;;------------------------------------------------------
;; NCL function for var utils 
;; Kai Zhang, MPI-M  
;;
;; list of funcitons: 
;; --- 
;;  function varsamedims(aa,bb)
;;  function varndims(aa)
;;  function varxthdim(aa,dimname)
;;  function varrvslat(aa)
;;  function varrvslat_name(aa,cc)
;;  function varstddim2d(aa)
;;  function xdimavg(aa,str)
;;  function vnmiss_1D(aa)
;;  function set_missing_1D(aa)
;;  function samedim_1D(aa,tm,tn)
;;  function var_pmel(aa,tm)
;;  function var_passend(aa,pid)
;;  function var_1DHH2D(aa,idebug)
;;  function var_2DHH1D(aa,idebug)
;;  function var_1H2D1215HH(aa,idebug)
;;------------------------------------------------------


;;******************************************************
;; to determine if two variables have the same
;; length for each dimension  
;;******************************************************
  undef("varsamedims")

  function varsamedims(aa,bb)

  begin

   da = dimsizes(aa)  
   db = dimsizes(bb)  

   if (da.eq.db) then 
      ;;print(" ") 
      ;;print("var --- dimension check done ")   
      ;;print(" ") 
   else
      print(" ") 
      print("var --- dimension check failed ")   
      print("exit") 
      exit
   end if 

   return(0)
 
  end



  undef("varndims")

;;------------------------------------------------------
;; to determine how many dementions a variable has  
;;------------------------------------------------------

  function varndims(aa)

  begin

   da = dimsizes(aa)  
   na = dimsizes(da)  

   return(na)
 
  end

;;------------------------------------------------------
;; based on the dimention name, determine which 
;; dimention it is  
;;------------------------------------------------------
  undef("varxthdim")
  function varxthdim(aa,dimname)
  begin
   na = varndims(aa) 
   id = -999 
   do ia = 0, na - 1 
      if(aa!ia.eq.dimname) then 
         id = na - ia 
      end if
   end do
   return(id)
  end

;;------------------------------------------------------
;; to reverse lat dimension
;;           input: aa - variable 
;;          output: is - status of processing   
;; function called: varndims and varxthdim  
;;------------------------------------------------------
  undef("varrvslat")
  function varrvslat(aa)
  begin

   na = varndims(aa) 
   id = varxthdim(aa,"lat") 
   is = -999 
 
   ;; 1,1 

   lyes = na.eq.1.and.id.eq.0

   if(lyes) then
      aa(      :) = aa(      ::-1)
      is = 0
   end if
 
   ;; 2,1 

   lyes = na.eq.2.and.id.eq.1 
 
   if(lyes) then 
      aa(    :,:) = aa(    ::-1,:) 
      is = 0 
   end if
   
   ;; 3,1 

   lyes = na.eq.3.and.id.eq.1 

   if(lyes) then 
      aa(  :,:,:) = aa(  :,::-1,:) 
      is = 0 
   end if

   ;; 4,1
   
   lyes = na.eq.4.and.id.eq.1 

   if(lyes) then 
      aa(:,:,:,:) = aa(:,:,::-1,:) 
      is = 0 
   end if

   ;; 3,2

   lyes = na.eq.3.and.id.eq.2 
   
   if(lyes) then 
      aa(  :,:,:) = aa(  ::-1,:,:) 
      is = 0 
   end if

   ;; 4,2

   lyes = na.eq.4.and.id.eq.2 
   
   if(lyes) then 
      aa(:,:,:,:) = aa(:,::-1,:,:) 
      is = 0 
   end if

   ;; 4,3
   
   lyes = na.eq.4.and.id.eq.3 

   if(lyes) then 
      aa(:,:,:,:) = aa(::-1,:,:,:) 
      is = 0 
   end if
   
   return(is)
  end

;;------------------------------------------------------
;; to reverse lat dimension
;;           input: aa - variable 
;;          output: is - status of processing   
;; function called: varndims and varxthdim  
;;------------------------------------------------------
  undef("varrvslat_name")
  function varrvslat_name(aa,cc)
  begin

   na = varndims(aa) 
   id = varxthdim(aa,"lat") 
   is = -999 
   
   ;;printVarSummary(aa)
   ;;print("na,id" + na + " " +id)
   
   if(na.eq.1.and.id.eq.1) then 
      aa(:) = aa(::-1) 
      is = 0 
   end if
   
   if(na.eq.2.and.id.eq.2) then 
      aa(:,:) = aa(::-1,:) 
      is = 0 
   end if
   
   if(na.eq.3.and.id.eq.2) then 
      aa(:,:,:) = aa(:,::-1,:) 
      is = 0 
   end if
   
   if(na.eq.4.and.id.eq.2) then 
      aa(:,:,:,:) = aa(:,:,::-1,:) 
      is = 0 
   end if
   
   if(na.eq.3.and.id.eq.3) then 
      aa(:,:,:) = aa(::-1,:,:) 
      is = 0 
   end if
   
   if(na.eq.4.and.id.eq.3) then 
      aa(:,:,:,:) = aa(:,::-1,:,:) 
      is = 0 
   end if
   
   if(na.eq.4.and.id.eq.4) then 
      aa(:,:,:,:) = aa(::-1,:,:,:) 
      is = 0 
   end if
   
   ;;if (is.eq.0) then 
   ;;print(cc + " is inverted in lat dim") 
   ;;end if
   
   return(is)
  end
  

;;------------------------------------------------------
;; under development  
;;------------------------------------------------------
  undef("varstddim2d")
  function varstddim2d(aa)
  begin
   aa!0 = "lat"
   aa!1 = "lon"  
   return(0)
  end

;;------------------------------------------------------
;; under development
;;------------------------------------------------------
  undef("xdimavg")
  function xdimavg(aa,str)
  begin

   na = varndims(aa)

   if (na.gt.5) then 
    print("dimsize larger than 5 !!! not supported !!!") 
    exit 
   end if 
  
   dimname = new(na,"string") 
   dimnwid = new(na,"integer") 

   ;; get dimension name and preset the new dimention index 

   do is = 0, na - 1 
    dimname(is) = aa!is
    dimnwid(is) = is  
    print(aa!is)     
   end do 

   lyes = any(dimname.eq.str) 

   if (.not.lyes) then 
    print("no "+str+" dimension !!!") 
    exit 
   end if 

   ;; get avg-dim index 

   id = varxthdim(aa,str)  

   ;; get the new dimention index 
   ;; increase index of dimension which is behind the avg-dim (id)  
   ;; eg. id=1, 4-1=3, 4-1=3 

   if (na.ge.2) then  
    do is = na - id, na - 1  
     dimnwid(is) = dimnwid(is) + 1  
    end do 
   end if 

   if (na.eq.1) then
    bb = dim_avg(aa)
   end if

   if (na.eq.2) then
    bb = dim_avg(aa($dimname(dimnwid(0))$|:,$str$|:))
   end if

   if (na.eq.3) then 
    bb = dim_avg(aa($dimname(dimnwid(0))$|:,$dimname(dimnwid(1))$|:,$str$|:))  
   end if  

   if (na.eq.4) then
    bb = dim_avg(aa($dimname(dimnwid(0))$|:,$dimname(dimnwid(1))$|:,$dimname(dimnwid(2))$|:,$str$|:))
   end if

   if (na.eq.5) then
    bb = dim_avg(aa($dimname(dimnwid(0))$|:,$dimname(dimnwid(1))$|:,$dimname(dimnwid(2))$|:,$dimname(dimnwid(3))$|:,$str$|:))
   end if

   return(bb)
  end


;;------------------------------------------------------
;; delete missing part of data 
;;------------------------------------------------------
  undef("vnmiss_1D")
  function vnmiss_1D(aa)
  begin

   idb = 1 

   ndim = dimsizes(dimsizes(aa))

   if (ndim.le.1) then

      a1 = aa

      if (idb.eq.0) then 
         print(a1@long_name+"")
      end if 

      a1 = where(a1.lt. -998.,-999,a1) 
      a1 = where(a1.gt.1.E+10,-999,a1) 


      a1@_FillValue = -999 

      dima1 = dimsizes(a1)

      lyes = all(ismissing(a1)) 

      if (lyes) then 
         print(" missing all values in: aa") 
         exit 
      end if 

      lyes = any(ismissing(a1)) 

      b1 = a1
      b1 = 0.  

      if (lyes) then

         is = 0

         do id = 0, dima1 - 1
            if (.not.ismissing(a1(id))) then
               b1(is) = a1(id)
               is = is + 1
            else
               ;;print("a1("+id+") = "+a1(id)+" is missing")  
            end if
         end do

         if (idb.eq.0) then 
            print("number of non-missing values: "+is)  
         end if  

         if (is.eq.0) then 
            print(a1@long_name+"")
            print(" a1 all missing values" )
            exit  
         end if 
         c1 = b1(0:is-1)

      else
         c1 = a1 
      end if

   else
      print(" vnmiss_1D needs 1D input !!!") 
      exit 
   end if

   delete(a1) 
   delete(b1) 

   return(c1)
  end


;;------------------------------------------------------
;; set missing values 
;; external: none
;;------------------------------------------------------
  undef("set_missing_1D")
  function set_missing_1D(aa)
  begin
   aa = where(aa.lt. -900.,-999,aa) 
   aa = where(aa.gt.1.E+10,-999,aa) 
   aa@_FillValue = -999 
   return(aa)
  end


;;------------------------------------------------------
;; delete missing part of data 
;;------------------------------------------------------
  undef("samedim_1D")
  function samedim_1D(aa,tm,tn)
  begin

   idb = 0 

   ndim = dimsizes(dimsizes(aa)) 

   tm = set_missing_1D(tm) 

   misstm = ismissing(tm) 

   nvam = dimsizes(tm) 
   nvan = dimsizes(tn)

   print(nvam+" "+nvan+" ")
   ;;print(misstm)
 
   if (ndim.le.1) then
 
      lyes = any(misstm) 

      cc = aa(0:nvan-1) 
      cc = 0.  
 
      if (lyes) then

         is = 0

         do id = 0, nvam - 1
            if (.not.misstm(id)) then
               cc(is) = aa(id)
               is = is + 1
            else
               ;;print("a1("+id+") = "+a1(id)+" is missing")  
            end if
         end do

         if (idb.eq.0) then 
            print("number of non-missing values: "+is)  
         end if  
 
         dd = cc(0:is-1)

      else
         dd = aa
      end if

   else
      print(" vnmiss_1D needs 1D input !!!") 
      exit 
   end if
 
   return(dd)
  end


;;------------------------------------------------------
;; PMEL data format
;; external: vnmiss_1D 
;;------------------------------------------------------
  undef("var_pmel")
  function var_pmel(aa,tm)
  begin

   idb = 0 

  ;tn = vnmiss_1D(tm)
   tn = tm

   bb = new((/dimsizes(tn(:))/),"float")

   bb!0 = "time"
   bb&time = tn
   bb@long_name = aa@long_name
   bb@units     = aa@units


   cc = samedim_1D(aa(:,0,0,0),tm,tn) 

   ;;cmd = prtsum(cc) 
   ;;cmd = prtsum(bb)  
 
   if (idb.eq.0) then  
      print(cc@long_name+"")
   end if 

   bb(0:dimsizes(cc)-1) = (/ cc /)

   return(bb)
  end



;;------------------------------------------------------
;; cut certain header lines of a var
;; and put them in the end
;;------------------------------------------------------
  undef("var_passend")
  function var_passend(aa,pid)
  begin

   nd = dimsizes(aa) 
   bb = new((/nd/),typeof(aa)) 
                                  ;; suppose nd  = 8 
   if (pid.ne.0) then             ;; suppose pid = 6  
     bb(0:nd-1-pid) = aa(pid:)    ;; bb(0:8-1-6) = aa(6:)  
    ;bb(nd-pid:)    = aa(0:pid-1) ;; bb(8-6)     = aa(0:6-1)  
     bb(nd-pid:)    = -999        ;; set to missing values 
   else
     bb = aa
   end if 
 
   return(bb) 
  end

;;------------------------------------------------------
;; convert 1D hourly data to 2D data 
;; which contains 24 elements in the 2nd dimension 
;;------------------------------------------------------
  undef("var_1DHH2D")
  function var_1DHH2D(aa,nh,idebug)
  begin

  nd = dimsizes(aa)
  nm = mod(nd,nh) 

  if(nm.ne.0) then 
     ns = (nd + nh - nm) / nh  
  else
     ns = nd / nh  
  end if

  if(idebug.eq.1) then 
     print("")
     print("--- var_1DHH2D ---")
     print("nd: "+nd) 
     print("nm: "+nm) 
     print("ns: "+ns) 
     print("")
  end if

  bb = new((/ns,nh/),typeof(aa))
  bb = -999 

  ix = 0 

  do is = 0,ns-1
  do it = 0,nh-1 
     if(ix.le.nd-1) then 
        bb(is,it) = (/ aa(ix) /) 
        ix = ix + 1 
        ;;if(idebug.eq.1) then 
        ;;   print("is: "+is+" it: "+it+" ix: "+ix)  
        ;;end if  
     end if
  end do 
  end do 

  return(bb)
  end


;;------------------------------------------------------
;; convert 2D hourly data 
;; which contains 24 elements in the 2nd dimension
;; to 1D data
;;------------------------------------------------------
  undef("var_2DHH1D")
  function var_2DHH1D(aa,idebug)
  begin

  nd = dimsizes(aa)
  ne = dimsizes(nd)
  ns = nd(0)
  nh = nd(1)
  nt = ns*nh

  if(idebug.eq.1) then

     if(ne.ne.2) then
        print("var_1H2D1215HH needs 2D var !!!")
        exit
     end if

     print("")
     print("--- var_2DHH1D ---")
     print("ne: "+ne)
     print("nd: "+nd)
     print("")

  end if


  bb = new((/nt/),typeof(aa))
  bb = -999

  ix = 0

  do is = 0,ns-1
  do it = 0,nh-1
     if(ix.le.nt-1) then
        bb(ix) = aa(is,it)
        ix = ix + 1
     else
        return 
     end if
  end do
  end do

  return(bb)
  end




;;------------------------------------------------------
;; select obs data from 12H to 15H 
;; for Cape Grim 
;; looking at afternoon mean radon concentrations 
;;------------------------------------------------------
  undef("var_1H2D1215HH")
  function var_1H2D1215HH(aa,idebug)
  begin

  nd = dimsizes(aa)
  ne = dimsizes(nd) 
 
  if(idebug.eq.1) then

     if(ne.ne.2) then 
        print("var_1H2D1215HH needs 2D var !!!")
        exit  
     end if 
   
     if(nd(1).ne.24) then 
        print("var_1H2D1215HH needs 24 elements in 2nd dimension !!!") 
        exit 
     end if 

     print("")
     print("--- var_1H2D1215HH ---") 
     print("ne: "+ne)
     print("nd: "+nd)
     print("")

  end if

  aa(:,:10) = -999 
  aa(:,15:) = -999 

  return(aa)
  end


;;------------------------------------------------------
;; select mod data from 12H to 15H 
;; for Cape Grim 
;; looking at afternoon mean radon concentrations 
;;------------------------------------------------------
  undef("var_3H2D1215HH")
  function var_3H2D1215HH(aa,idebug)
  begin

  nd = dimsizes(aa)
  ne = dimsizes(nd) 
 
  if(idebug.eq.1) then

     if(ne.ne.2) then 
        print("var_3H2D1215HH needs 2D var !!!")
        exit  
     end if  

     print("")
     print("--- var_3H2D1215HH ---") 
     print("ne: "+ne)
     print("nd: "+nd)
     print("")

  end if

  aa(:,:3) = -999 
  aa(:,6:) = -999 

  return(aa)
  end

;;------------------------------------------------------
;; select obs data from 12H to 15H
;;------------------------------------------------------
  undef("var_1HH1215")
  function var_1HH1215(aa,idebug)
  begin

   a1 = var_1DHH2D(aa,24,1)   ;; 1D to 2D
   a2 = var_1H2D1215HH(a1,1)  ;; select 12-15H local time
   a3 = var_2DHH1D(a2,1)      ;; 2D to 1D

  return(a3)
  end



;;------------------------------------------------------
;; select mod data (3hourly) from 12H to 15H
;;------------------------------------------------------
  undef("var_3HH1215")
  function var_3HH1215(aa,idebug)
  begin

   a1 = var_1DHH2D(aa,8,1)    ;; 1D to 2D
   a2 = var_3H2D1215HH(a1,1)  ;; select 12-15H local time
   a3 = var_2DHH1D(a2,1)      ;; 2D to 1D

  return(a3)
  end


;;------------------------------------------------------
;; P2 percent  
;;------------------------------------------------------
  undef("var_p2")
  function var_p2(aa,bb,idebug)
  begin

  ic = 0
  ia = 0

  do i = 0, dimsizes(aa)-1
     if(.not.ismissing(aa(i)).and..not.ismissing(bb(i))) then
        if(aa(i) .ge. 0.00001) then
           if(bb(i)/aa(i).le. 0.5 .or. bb(i)/aa(i).ge. 2.0 ) then
               ic = ic + 1
           end if
           ia = ia + 1
        end if
     end if
  end do

  if(ia.ne.0) then
    fy = 100. - int2flt(ic) / int2flt(ia) * 100.
  else
    fy = 0.
  end if

  return(fy)
  end

;;------------------------------------------------------
;; mask negative  
;;------------------------------------------------------
  undef("var_maskneg")
  function var_maskneg(aa,idebug)
  begin

  aa = where(aa.lt.0.,-999.,aa)

  return(aa)
  end


;;------------------------------------------------------
;; mask for DJF  
;;------------------------------------------------------
  undef("var_onlyDJF")
  function var_onlyDJF(aa,idebug)
  begin

  ndm = varndims(aa) 
  idt = varxthdim(aa,"time") 

  bb = aa 

  ist = 0
 
  if(ndm.eq.2) then 
     if(idt.eq.1) then 
        bb(:,2:10) = -999 ;; only DJF
     end if 
     if(idt.eq.2) then 
        bb(2:10,:) = -999 ;; only DJF
     end if 
     ist = 1 
  end if

  if(ndm.eq.1) then
     if(idt.eq.1) then
        bb(2:10) = -999 ;; only DJF
     end if
     ist = 1
  end if

  if(ist.eq.0) then 
     print("var_onlyDJF only support 1D,2D var") 
  end if 

  return(bb)
  end


;;------------------------------------------------------
;; mask for JJA
;;------------------------------------------------------
  undef("var_onlyJJA")
  function var_onlyJJA(aa,idebug)
  begin

  ndm = varndims(aa)
  idt = varxthdim(aa,"time")

  bb = aa

  ist = 0 

  if(ndm.eq.2) then
     if(idt.eq.1) then
        bb(:,1:4) = -999 ;; only JJA
        bb(:,8: ) = -999
     end if
     if(idt.eq.2) then
        bb(1:4,:) = -999 ;; only JJA
        bb(8:, :) = -999
     end if
     ist = 1 
  end if

  if(ndm.eq.1) then
     if(idt.eq.1) then
        bb(0:4) = -999 ;; only JJA
        bb(8: ) = -999
     end if
     ist = 1
  end if

  if(ist.eq.0) then
     print("var_onlyJJA only support 1D,2D var")
  end if

  return(bb)
  end


;;------------------------------------------------------
;; mask for MAM
;;------------------------------------------------------
  undef("var_onlyMAM")
  function var_onlyMAM(aa,idebug)
  begin

  ndm = varndims(aa)
  idt = varxthdim(aa,"time")

  bb = aa

  ist = 0 

  if(ndm.eq.2) then
     if(idt.eq.1) then
        bb(:,0:1) = -999 ;; only MAM
        bb(:,5: ) = -999
     end if
     if(idt.eq.2) then
        bb(0:1,:) = -999 ;; only MAM
        bb(5: ,:) = -999
     end if
     ist = 1 
  end if 

  if(ndm.eq.1) then
     if(idt.eq.1) then
        bb(0:1) = -999 ;; only MAM
        bb(5: ) = -999
     end if
     ist = 1
  end if

  if(ist.eq.0) then
     print("var_onlyMAM only support 1D,2D var")
  end if

  return(bb)
  end


;;------------------------------------------------------
;; mask for SON 
;;------------------------------------------------------
  undef("var_onlySON")
  function var_onlySON(aa,idebug)
  begin

  ndm = varndims(aa)
  idt = varxthdim(aa,"time")

  bb = aa
  ist = 0 

  if(ndm.eq.2) then
     if(idt.eq.1) then
        bb(:,0:7) = -999 ;; only SON
        bb(:,11:) = -999
     end if
     if(idt.eq.2) then
        bb(0:7,:) = -999 ;; only SON
        bb(11:,:) = -999
     end if
     ist = 1 
  end if 

  if(ndm.eq.1) then
     if(idt.eq.1) then
        bb(0:7) = -999 ;; only SON
        bb(11:) = -999
     end if
     ist = 1
  end if

  if(ist.eq.0) then
     print("var_onlySON only support 1D,2D var")
  end if

  return(bb)
  end

;;------------------------------------------------------
;; min in fortran
;;------------------------------------------------------
  undef("formin")
  function formin(a,b)
  begin

  if (a.ge.b) then
     c = b
  else
     c = a  
  end if 

  return(c)
  end

;;------------------------------------------------------
;; max in fortran
;;------------------------------------------------------
  undef("formax")
  function formax(a,b)
  begin

  if (a.le.b) then
     c = b
  else
     c = a
  end if

  return(c)
  end









